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Abstract 

A  least-squares  technique  is  developed  for  identifying  unknown  parameters  in  a 
coupled  system  of  nonlinear  size-structured  populations.  Convergence  results  for  the 
parameter  estimation  technique  are  established.  Ample  numerical  simulations  and 
statistical  evidence  are  provided  to  demonstrate  the  feasibility  of  this  approach. 
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1  Introduction 


A  typical  direct  problem  for  structured  populations  is  to  use  the  knowledge  of  underlying 
mechanism  at  individual  level  such  as  growth,  mortality  and  reproduction  rates  to  deduce 
the  behavior  at  population  level.  This  approach  has  been  extensively  studied  for  many  kinds 
of  models  which  include  structured  and  non-structured  populations.  In  practice,  however, 
our  knowledge  of  the  vital  rates  may  be  incomplete  [40].  In  fact,  in  many  animal  and  plant 
populations  the  processes  at  the  individual  level  are  not  accessible  to  direct  observation 
[47].  For  example,  for  nonlinear  structured  models  the  dependence  of  reproduction  and 
mortality  rates  on  the  total  population  is  sometimes  completely  unknown  [37] .  Even  for  linear 
structured  models,  one  may  not  be  able  to  obtain  the  exact  dependence  of  the  vital  rates 
on  the  age  or  size  structure  [40].  In  these  cases,  one  resorts  to  an  inverse  problem  approach, 
namely  to  use  knowledge  about  the  behavior  at  the  population  level  (e.g,  observations  of 
total  population  numbers)  to  deduce  the  underlying  mechanisms  at  the  individual  level. 

In  recent  years  many  researchers  have  focused  their  attention  on  developing  methodologies 
for  solving  inverse  problems  governed  by  structured  population  models  (e.g,  [l]-[3],  [12]- [17] , 
[19]- [23],  [25]- [34],  [40]- [49]).  In  what  follows,  we  briefly  review  some  of  the  recent  work  on 
such  inverse  problems.  For  age-structured  population  models,  several  approaches  have  been 
developed  to  recover  unknown  individual  vital  rates.  For  example,  in  [40,  43]  a  fixed  point 
iterative  technique  was  developed  to  determine  the  death  rate  from  census  data  on  the  age 
distribution  of  the  population.  Therein,  conditions  on  the  data  are  given  that  lead  to  a 
unique  solution.  In  [26]  the  authors  formulated  the  inverse  problem  as  an  operator  equation 
and  the  least  squares  method  is  then  used  to  compute  its  solution.  Due  to  the  ill-posedness 
of  the  problem,  a  regularization  technique  was  considered.  In  addition,  the  authors  prove 
that  the  resulting  scheme  has  a  convergence  rate  of  Holder  type.  However,  no  numerical 
results  were  reported.  A  least  squares  approach  was  also  adopted  in  [19]  for  a  nonlinear  age- 
structured  population  model  to  estimate  unknown  coefficients  from  a  set  of  fully  discrete 
observations  of  the  population.  Although  the  convergence  of  the  computed  minimizers  to  a 
minimizer  of  the  least  squares  problem  was  established  and  numerical  results  were  presented, 
for  many  real  populations  it  is  generally  difficult  to  obtain  discrete  observations  with  respect 
to  age,  whereas  other  quantities  such  as  total  population  number  are  easily  obtained.  In  [25] 
a  model  describing  the  evolution  in  time  of  size/age  structured  population  was  considered. 
A  moving  finite  element  method  was  used  to  study  the  identification  problem  for  such  a 
model.  Convergence  results  for  the  parameter  estimation  technique  were  reported.  In  [30], 
by  writing  a  linear  age-structured  model  using  the  cumulative  formulation  approach  (see  e.g., 
[24]),  the  authors  studied  the  inverse  problem  of  identifying  the  birth  and  death  rates  from 
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data  on  the  total  population  size  and  the  cumulative  number  of  births.  They  also  provided 
conditions  on  the  data  that  guarantee  the  uniqueness  of  the  solution  to  the  inverse  problem. 

For  size-structured  population  models,  the  least  squares  approach  has  been  often  used  for 
parameter  identification.  For  example,  it  was  used  in  [15,  16]  to  estimate  the  growth  rate  dis¬ 
tribution  in  a  linear  size-structured  population  model.  A  similar  technique  was  subsequently 
applied  to  a  semilinear  size-structured  model  in  [34]  where  the  mortality  rate  depends  on 
the  total  population  due  to  competition.  In  [2]  an  inverse  problem  governed  by  a  phyto¬ 
plankton  aggregation  model  was  studied.  Convergence  and  numerical  results  for  identifying 
the  coagulation  kernel  were  provided.  Later,  this  technique  was  extended  to  identify  pa¬ 
rameters  in  a  size-structured  population  model  in  [1,  3]  where  all  the  individual  vital  rates 
(growth,  mortality  and  reproduction)  depend  on  the  total  population  level.  Therein,  these 
parameters  are  identified  from  a  set  of  observations  corresponding  to  the  total  population 
number.  A  finite  difference  method  was  then  used  to  approximate  the  infinite  dimensional 
problem.  Convergence  results  for  the  computed  parameter  estimates  to  the  true  parameter 
were  established.  To  our  knowledge,  [3]  was  the  first  paper  to  provide  convergence  results 
for  parameter  estimates  when  the  growth  rate  is  a  nonlinear  function  of  the  total  population 
(i.e.,  the  size-structured  model  is  represented  by  a  quasilinear  first  order  hyperbolic  initial 
boundary  value  problem). 

In  this  paper  we  extend  the  discussion  in  [3]  to  the  following  coupled  system  of  quasilinear 
size-structured  populations  model: 

u{  +  (g\x,  Pit ;  q))uI)x  +  m/(x,  P(t;  q))u:  =  0,  ( x ,  t)  E  (0,  L\  x  (0,  T], 

N  rL 

gI(Q,P(t;q))uI(0,t;q)  =  CI(t)  +  '^2  7 I,J pJ (x,  P(t;  q))uJ (x,t;  q)dx,  le(0,T],  (1.1) 

j=iV° 

uI(x10]q)  =  uI'0(x),  iG  [0,1]. 

Here  q  =  {q1,  q2, . . . ,  qN )  with  q1  =  (g1 ,  m1 ,  (31 ,  C1),  I  =  1,  2, . . . ,  N,  the  parameters  to  be 
identified.  The  function  uJ(x,  t;q),  I  =  1,  2, . . . ,  N,  is  the  parameter-dependent  size  density 
(number  per  unit  size)  of  individuals  in  the  Jth  population  having  size  x  at  time  t,  and 

N  .L 

=  uJ(x,t;q)dx  (1.2) 

j=i  Jo 

is  the  total  population  at  time  t.  The  function  g1  denotes  the  growth  rate  of  an  individual 
in  the  Jth  population,  m1  denotes  the  mortality  rate  of  an  individual  in  the  Jth  population, 
and  p1  is  the  reproduction  rate  of  an  individual  in  the  Jth  population.  The  function  C1 
represents  the  inflow  rate  of  the  Jth  population  of  zero-size  individuals  from  an  external 
source  (e.g.,  in  a  tree  population  model  seeds  moved  by  wind). 
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The  model  (1.1),  which  was  developed  by  the  authors  in  [4],  is  a  generalization  of  several 
size-structured  population  models  (usually  referred  to  as  structured  models  with  rate  dis¬ 
tributions)  which  have  been  investigated  in  [14,  15,  16,  34],  Motivated  by  the  fact  that,  in 
addition  to  observable  characteristics  such  as  age  or  size  of  the  individuals,  non-observable 
genetic  characteristics  may  often  play  a  crucial  role  in  the  development  of  the  individuals, 
researchers  in  [14]  presented  the  first  such  generalization  of  the  classical  Sinko-Streifer  model. 
This  model,  which  is  a  linear  version  of  (1.1),  has  vital  individual  rates  that  are  independent 
of  the  total  population  and  distributed  over  an  an  infinite-dimensional  admissible  parame¬ 
ter  space  with  a  probability  measure.  It  was  shown  through  numerical  simulations  in  [14] 
that  there  is  a  crucial  difference  between  the  dynamics  of  distributed  rate  size-structured 
population  models  and  the  classical  Sinko-Streifer  models.  In  particular,  the  classical  Sinko- 
Streifer  model  cannot  have  dispersion  of  the  density  of  the  population  in  age  or  size  except 
under  biologically  unreasonable  conditions  on  the  growth  rate  [15].  That  is  why  the  classi¬ 
cal  Sinko-Streifer  models  are  in  conflict  with  field  data  collected  by  experimental  biologists. 
These  data  sets  show  that  a  population  with  unimodal  distribution  evolves  into  a  bimodal 
distribution  (see  [14]  and  [41]).  In  [17]  the  authors  used  least  squares  approach  to  fit  these 
distributed  rate  models  to  data  obtained  in  [14],  The  resulting  good  fit  indicates  that  the 
need  for  such  modification  is  crucial  if  these  models  were  to  be  used  as  prediction  tools. 

In  addition  to  extending  the  theory  in  [3]  to  the  coupled  quasilinear  system  (1.1),  a  main 
novelty  of  our  current  research  is  that  we  report  on  extensive  numerical  simulations.  These 
simulations  are  then  used  to  obtain  statistical  results  (in  the  form  of  confidence  intervals) 
which  provide  solid  evidence  on  the  feasibility  of  this  approach.  It  is  worth  pointing  out 
that  with  the  exception  of  [28]  the  above-mentioned  articles  do  not  report  on  any  statistical 
studies. 

As  the  use  of  numerical  methods  for  estimating  functional  parameters  becomes  more 
widely  accepted  in  the  biological  sciences,  it  is  becoming  increasingly  important  for  inves¬ 
tigators  to  support  the  efficacy  of  proposed  numerical  algorithms  with  not  only  numerical 
simulation  results  but  also  confidence  intervals  on  estimated  parameters.  This  can  be  done 
by  calculating  standard  errors  in  a  number  of  sophisticated  ways  (e.g.,  pointwise  confidence 
intervals  or  bands  as  in  [38,  39,  48],  uniform  bands  [32],  simultaneous  confidence  bands  [31], 
etc.).  Here  we  simply  compute  the  pointwise  standard  errors  using  the  pointwise  sample 
variances  from  a  large  (1000)  number  of  inverse  problem  simulations.  While  in  our  efforts  we 
emphasize  (regularized)  ordinary  least  square  estimators,  the  ideas  and  methods  presented  in 
this  paper  can  readily  be  used  with  maximum  likelihood  estimators  as  well  as  other  standard 
estimators  found  in  the  statistical  literature. 

It  is  also  worth  noting  another  connection  between  statistical  methods  and  our  efforts 
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in  this  paper.  The  models  we  use  here  involve  a  form  of  “mixing”  distributions  found  in 
the  literature  on  mixed  effects,  random  effects  or  hierarchial  methods  (see  for  example, 
[20,  21,  22,  35,  36,  46]).  However  the  models  we  investigate  entail  mixing  that  cannot  be 
decoupled  into  individual  dynamics  and  thus  result  in  fully  coupled  dynamics  (see  our  closing 
remarks  in  Section  4). 

By  a  weak  solution  to  problem  (1.1)  we  mean  a  bounded  and  measurable  function 
u{x,  t;  q )  =  (u1(a:,  t;  q ),  u2{x,  £;  q), . . . ,  uN(x ,  £;  q))  satisfying 

r*L  r>L 

/  u1  (x,t;  q)cp(x,t)dx  —  /  u\x,  0;  q)ip{x,  0)dx 

Jo  Jo 


(u1  ips  +  gV <px  —  m1uI(p)dx  ds 


o  Jo 


(1.3) 


+  [  ^(o,  s)  ( c7(s) + y 

Jo  \  J=1  Jo 


7 1,J /3J (x,  P(s-,  q))uJ (x,  S]  q)dx  )  ds 


for  t  G  [0,  T],  I  =  1,  2, . . . ,  N,  and  every  test  function  ip  e  C1([0,  L]  x  [0,  T]). 
We  first  impose  a  condition  on  the  initial  data:  for  any  /  =  1,2, ...  ,N 


(HI)  uIfi  G  BV[0,L\  and  uIfi(x)  >  0. 

N 

Then  let  B  =  B1  with  B1  =  C1([0,  L];  Cb[0,  oo))  x  C&(0)  x  Cb(fi)  x  C[0,T\,  where 

i=i 

fl  =  [0,  L\  x  [0,  cx))  and  C&( fl)  denotes  the  space  of  uniformly  bounded  continuous  functions 

on  O.  We  assume  that  our  admissible  parameter  space  Q1  is  a  compact  subset  of  B1  satisfying 

(H2)-(H5)  below. 

(H2)  p!(x,  P)  is  a  nonnegative  Lipschitz  continuous  function  in  x  and  P  with  a  Lipschitz 
constant  Li.  Furthermore,  ^{x,  P)  <  where  uj\  is  a  positive  constant. 

(H3)  m7(a:,P)  is  a  nonnegative  Lipschitz  continuous  function  in  x  and  P  with  a  Lipschitz 
constant  L2.  Furthermore,  m7(a;,P)  <  uj2,  where  u2  is  a  positive  constant. 

(H4)  g1  (x,  P )  is  twice  continuously  differentiable  with  respect  to  x  and  satisfies  | g\x,  P)|  + 
| gi(x,  P )  |  +  \gTxx(x,  P )  |  <  cu3,  where  u3  is  a  positive  constant.  Furthermore,  g\x,  P)  >  0 
for  x  E  [0,  L)  and  g\L,P )  =  0,  and  g\x,P)  and  gl.{x,P)  are  Lipschitz  continuous  in 
P  with  a  Lipschitz  constant  L3. 


(H5)  C1  (t)  is  a  nonnegative  Lipschitz  continuous  function  with  a  Lipschitz  constant  L4 
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N 


Let  Q  =  Q1 ,  then  Q  is  a  compact  subset  of  B. 


1=1 

Depending  on  the  values  of  the  constants  0  <  y7’J  <  1,  the  model  (1.1)  may  have 
two  different  interpretations.  If  7 1,1  =  1  and  y7,7  =  0,  J  ^  ./,  the  model  represents  the 
dynamics  of  several  populations  competing  for  common  resources.  On  the  other  hand,  if 
y7,7  >  0,  J,  J  —  1, 2, . . . ,  N,  then  the  model  may  describe  the  dynamics  of  one  population 
consisting  of  N  subpopulations,  each  with  its  own  characteristics.  Hence,  y7,7  represents 
the  probability  that  an  individual  of  the  Jth  subpopulation  will  reproduce  an  individual  of 
the  Jth  subpopulation.  Therefore,  two  different  ways  for  observing  data  will  be  considered. 
These  lead  to  the  following  two  different  least-squares  functionals  to  be  minimized:  The 
first  one  is  based  on  the  assumption  that  the  model  (1.1)  describes  N  different  competing 
populations.  Hence  observations  ZIk  which  correspond  to  the  total  number  of  individuals  in 
the  Jth  population  at  time  4  are  assumed  to  be  available  (this  case  corresponds  to  y7,7  =  1 
and  y7,7  =  0,  J  ^  J).  We  define  the  least-squares  cost  functional  for  this  case  to  be 

,2 


J{q)  = 


1  k 


log 


u\x,  4;  q)dx  +11-  log(ZI>k  +  1) 


(1.4) 


which  is  minimized  over  Q.  The  other  case  assumes  that  (1.1)  models  one  species  which 
has  been  divided  into  N  not  readily  distinguishable  subpopulations.  In  this  case,  we  assume 
that  we  can  only  observe  aggregate  data  Zkl  the  total  number  of  individuals  at  time  tk  (this 
case  corresponds  to  y7,J  >  0,  J,  J  —  1,  2, . . . ,  N).  We  define  the  least-squares  cost  functional 


J(q)  = 


y. 

k 


u^x,  4;  q)dx  +  1  )  -  log (Zk  +  1) 


(1.5) 


which  is  minimized  over  Q. 

We  remark  that  minimizing  (1.4)  over  Q  is  equivalent  to  the  maximum  likelihood  esti¬ 
mation  of  q  if 

ei,k  =  log  (^J  u\x,  4;  q)dx  +  -  log  (ZIjk  +  1) 

are  i.i.d.  normal,  and  minimizing  (1.5)  over  Q  is  equivalent  to  the  maximum  likelihood 
estimation  of  q  if 


Cfc  =  log 


uI(x,tk;q)dx+  1  -  log(Zfc  +  1) 


are  i.i.d.  normal. 

The  paper  is  organized  as  follows.  In  Section  2,  we  present  a  finite  difference  scheme 
for  computing  the  solution  of  (1.1)  and  then  provide  convergence  results  for  the  parameter 
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estimation  technique.  In  Section  3,  we  give  ample  numerical  and  statistical  results.  Some 
concluding  remarks  are  made  in  Section  4. 


2  Approximation  Scheme  and  Convergence  Theory 

The  following  notation  will  be  used  throughout  the  paper:  Ax  =  L/n  and  A t  —  T/l  denote 
the  spatial  and  time  mesh  size,  respectively.  The  mesh  points  are  given  by  Xj  =  j  Ax,  j  = 
0, 1,  2, . . . ,  n  and  tk  =  kAt,  k  =  0, 1, 2 ,...,/.  We  denote  by  u*(q)  and  Pk(q )  the  finite 
difference  approximation  of  u1  (xj,tk,  q)  and  P(tk,q),  respectively,  and  we  let 

g]'k  =  gI(xj,Pk(q)),  p!’ k  =  pI(xj,Pk(q)), 


m 


**  =  m\xj,Pk(q)),  and  C 1)k  =  Cy(4). 


I,k 


We  define  the  difference  operator 


D~Mk)  = 


I,k  I,k 

ui  -  Vi 


Ax 


1  <j  <  n 


and  the  i1,  £°°  and  the  B V  norms  of  u1,K  by 


I,k 


\uI,k\\i  =  ^2  \uj,k\^x’ 

3= 1 


i  T  h  ii  i  T  h  \  II  T  h  n 

\u  '  ho  =  max  «■’  ,  \\u  '  \\bv  = 

0 <j<n  J 


X  X,T(«f  )|Ax. 

3= 1 


We  then  discretize  the  partial  differential  equation  in  (1.1)  using  the  following  implicit  finite 
difference  approximation 

I,k+1 (  \  7,/c/  \  I.k  7.&+1/  \  I.k  I,k-\-l/  \ 

(?)-«,■  W  ,  9j  ui  W**9i- iui-i  W) 


At 


+ 


Ax 


+  mI*uIj’k+1{q)  =0,  1  <  j  <  n, 


N  n 


d*<-k+I  w  =  c* +Y.Y.  <«) 

N  n 

Pk+\q)  =  J2J2uI3,k+1^Ax 

i=i  j= i 


(2.1) 


with  the  initial  condition 


i  rjA  x 

uIj’°  =  —  uIfi(x)dx,  j  =  1,  2, . . . ,  n. 

Ax  i)Ax 


If  we  define 


d1/  =  1  +  ^ g)*  +  Atm**  j  =  1,  2, . . . ,  n,  I  =  1,  2, . . . ,  A, 
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then  (2.1)  can  be  equivalently  written  as  the  following  system  of  linear  equations  for 


uk+1  (q)  =  Uq’k+1  (q) ,  < 'k+l  (q),...,  u]f+1  (q) ,  (q) ,  (q) , , , . ,  u2’k+1  (q),..., 


i,fe+i/  \  i,fc+i 


l,fc+l/  \  2,fc+l  / 


2,fe+l/ 


2,fe+l  / 


uZ’k+1(q),U?’k+\q),...,U^+\q) 


N,k+1 1 


Akuk+\q)  =  f\q). 


t>Nx(n+l) 


where 


/*(«)  = 


N  n 

cl*  +  EE  T'-yf  «f  (?)  Ai,  «!•»(«), ....  nj*(»), 

J= 1  i=i 
AT  n 

C'2’*  +  7  2’Jp-’kufk(q)Ax,  u2{k(q), u2/(q), 

J=  l  l=i 

/T1 

AT  n 


C 

»'‘+EEW(!)  a*, 

N,kt  \ 

A  (q), 

J= 

=1 1=1 

and  Ak  is  the  following  block  diagonal  matrix: 

(  Ahk  0  0  •  •  • 

°  \ 

o 

ci 

O 

0 

Ak  = 

0  0  A3’k  ■  ■  ■ 

0 

\  0  0  0  ••• 

^  . 

j?  ; 

?r 

with  the  lower  triangular  matrix 

( 

I,k 

% 

0  0  ••• 

0 

At  I,k  d,k  n 

—xs<,  <*.  0  ■■■ 

0 

AI,k  = 

0 

At  I,k  jl,k 

Axg'  " 

0 

At  j 

V 

0 

0  0 

Ax9n 

0  \ 
0 

0 


<’fc 


y 


(2.2) 


Note  that  using  the  assumptions  on  our  parameters  one  can  easily  show  that  equation 
(2.2)  has  a  unique  solution  satisfying  uk+1(q)  >  0,  k  —  0, 1, . . .  7 1  —  1. 

The  above  approximation  can  be  extended  to  a  family  of  functions  {U^x  At(x,t;q)}  de¬ 
fined  by 

UL  i((x,i;«)  =  «'•*(«)  for  (x,t)  £  x  [«*_!,«*), 

(2.3) 

j  =  1,  2, . . . ,  n,  k  —  1, 2, . . . ,  1,  I  =  1,  2, . . . ,  N. 


Since  our  parameter  set  is  infinite  dimensional,  a  finite  dimensional  approximation  of  the 
parameter  space  is  also  necessary  for  computing  minimizers.  To  this  end,  we  consider  the 
following  finite-dimensional  approximations  of  (1.4)  and  (1.5),  respectively: 


Jax,aM)  =  log  (  [ 

/  /.  \Jo 


UAxAt(X^^  (-l)dx  +  1  -  log (ZI,k  +  1) 


(2.4) 


and 

(2.5) 

each  of  which  is  minimized  over  Qm,  a  compact  finite-dimensional  approximation  of  the 
parameter  space  Q.  In  order  to  establish  the  convergence  results  for  the  parameter  estimation 
technique,  we  use  a  similar  approach  to  that  in  [3],  which  is  based  on  the  abstract  theory  in 
[18]- 


Ax,A ti^O)  ^  ^ 


log  El  UAxAt(X’ tk ;  ?) dx  +  1  )  -  log (zk  +  1) 


Theorem  2.1  Let  ql  =  (ql'l,q2’1 
Axi}  A ti  — »  0  as  i  — >  oo.  Let 


qN'1)  and  suppose  that  for  each  I ,  q1’1  — >  q1  in  Q 1  and 


I,i 


I  A.r  .a/.  t- ql)  =  (U^XiAt.(x,  t;  ql),  U'£XiAt.(x,  t;  q1),...,  U%XiAt.(x,  t;  q1)) 


denote  the  solution  of  the  finite  difference  scheme,  and  let 

u(x,  t;  q)  =  (u1^,  t; q),u2(x ,  t;q),...,  uN(x,  t;  q)) 
be  the  unique  weak  solution  of  our  problem  with  initial  condition 

u°(x)  =  (ul,0(x),u2’°(x), . . . ,  uN,0(x)) 

and  parameter  q,  then  At.(x,t;  ql)  — >  u1  (x,t;q)  mT1(0,  L)  uniformly  in  t  e  [0,T], 
Proof.  Dehne  ufk'1  =  ufk{ql).  From  the  fact  that  Q1  is  compact  and  the  results  of  [4],  there 

N 

exist  positive  constants  c i,  c2,  c3,  C4  such  that  for  each  I  =  1,  2, . . . ,  N,  we  have  £  \\uIAi\\i  < 


1=1 


Cl,  ||u/,fc’l||oo  <  C2,  II uIXl\\BV  <  c3  and  ^ 

3= 1 


u 


I,r,i 


U 


I,s  ,i 


At i 


A Xi  <  c4(r  —  s),  where  r  >  s.  Thus, 


for  each  /  there  exists  a  BV([0,L\  x  [0,  T])  function  uJ(x,t )  such  that  U^x.  Af.  (a;,  t;  ql)  — > 
in  £1(0,L)  uniformly  in  t.  Hence,  from  the  uniqueness  of  bounded  variation  weak 
solutions  stated  in  [4],  we  only  need  to  show  that  u(x,t )  =  (u1(x,t),u2(x,t), . . .  ,uN(x,t)) 
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is  the  weak  solution  corresponding  to  the  parameter  q.  To  this  end,  we  multiply  the  first 
equation  of  (2.1)  by  (pk+1  =  (p(xj,tk+ 1),  where  <p  G  C1([0,L]  x  [0,  T]),  to  obtain 

I,k+l,i  fc+1  I,k,i  k  jfc+1  k  I.k.i  I,k+l,i  fc+1  I,k,i  I,k+l,i  k+l 

u I  fj_  ZJh  n  _  +  9i  Vi  Z  gJA  MM 

A  ti  J  A  ti  Axi 

ink+l  —  </A+1 

I,k-\-l,i  r  j  r  j—1  I,k,i  I,k+l,i  fc+1 

-Sj -1  ",  :  '  at  +  ",  ^  =  °- 


Ax* 


Multiplying  the  above  equality  both  sides  by  Ax,  At,  and  summing  over  j  =  l,2,...,n, 
k  =  0, 1, ...,/  —  1,  we  find 


n  /-I  n  fc+1  _  ,  Jfc 

(  I,l,i  l  ho ,i  o\  A  hM  V j  Vj 

2^[ui  Vj-Uj  <Pj)  Axi  —  2_^ 


Ua 


,'=i 


k= 0  j=l 


1  1  nI,k,i  I, k+l, i  k+l  _  I,k,i  I, k+l, i  k+l 
iin  un  Yn  i)  0  u0  YO 


E 

k= 0 
l—l  n 


Ax; 


At; 


Ax,  At, 


l—l  n 


AxjAt; 


fc+i _  fc+i 

-EE  9j  I  I  ''  '  A/r.  J  A-AA/,:*  EE<M‘+1'  Vj+1AxjAfi  =  0. 


fc=o  j=i  fc=o  j=i 

Since  grnk,t  =  0  and  qI,r  — >  q1  as  i  — »  oo  in  Q7,  passing  to  the  limit  we  have 


u\x,  t)ip(x,  t)dx  —  /  u7(x,  0)g?(x,  0)<Tr 


0  ^0 


(w7y7  +  g7w7 \px  —  m7w7 \p)  dx  ds 

N  ,L 


+  f+  :m  fc'w+f:/' 

,0  V  ,0 


7 1,J /3J (x,  P(s))uJ (x,  s)dx  ds. 


Thus,  w(x,  f)  is  the  weak  solution  corresponding  to  the  parameter  g.  ■ 

Since  the  logarithm  function  is  continuous  on  [l,oo),  as  an  immediate  consequence  of 
Theorem  2.1,  we  obtain  the  following: 

Corollary  2.2  Let  UAx,At  denote  the  numerical  solution  of  (2.1)  with  parameter  ql  — »•  q  and 
Axi,  Ati  — ■>  0.  Then 

JAxiAuW)  -»•  J{q),  as  i  — >  oo. 


In  the  next  theorem,  we  establish  the  continuity  of  the  approximate  cost  functional,  so 
that  the  computational  problem  of  finding  approximate  minimizer  is  well-posed. 

Theorem  2.3  Let  Ax  and  At  be  fixed.  For  each  q1  G  Q1 ,  let  UAxAt(x,t;q)  denote  the 
solution  of  the  finite  difference  scheme,  and  qI)l  — >  g7  as  i  — >  oo  in  Q1 ,  then  UAx  At(x,  t;  ql )  — > 
UAx  At(x,  g)  as  i  ^  oo  in  £1(0 ,L)  uniformly  in  t  G  [0,T], 
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Proof.  Define  {uk’  ’*}  and  {uk’  }  to  be  the  solution  of  the  finite  difference  scheme  with 
parameter  q1  and  q,  respectively.  Let  =  Uj’k,t  —  Uj,k ,  then  v j,k'1  satisfies  the  following: 


I.k.i 

Va'  ~  Va 

- - x - —  +  D 

At 


h 


Pk^k+U  -  gI(xj,  Pk)uIM1 


+mI’i(xj,Pk’i)vIj’k+1’i  +  -  m^Xj^pb)]  u]'k+1  =  0, 


for  1  <  j  <  n,  and 


gJ,*(0  Pk,i)uI0’k+1,i  -  " Ifn  E>fc'V'Lfc+1 


V(o,PfeK 

N  n 


=  C^(tk)  -  C\tk )  +  Y,  E  P^vf^Ax 

J=  1  3= 1 


JV  n 


(2.6) 


(2.7) 


+  EE 7'''  !**)  ~  ?(*»> P‘)]  “f 

J=1  J  =  1 

where  PA:,?  denotes  Pk{ql).  Multiplying  both  sides  of  (2.6)  by  sgn(nj’fc+1’*)Ax  and  summing 
over  j  —  1,  2, . . . ,  n,  we  obtain 

II  11^7,  k,  ill 


At 


<  ~^2Dh  9I,l(xj,Pk't)Uj’k+1*  -  g\xjlPk)uIj'k+1  sgn(Vj,k+1,i) Ax 

3  = 1 


(2.8) 


-J2ml,i(xj,pk,i) 

3= 1 


j 


Ax 


E  W'Kx^P**)  -m7(^,Pfc)]  ^’fe+1sgn(vj’fe+M)Ax. 
1=1 

Using  the  fact  for  any  a,j  with  a,j  >  0,  j  =  0, 1,  2, . . . ,  n,  we  have 

n 

E  (aA)sgn(&j)Aa;  >  an|6n|  -  a0|6o|, 
j=i 
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we  obtain 


~YDh  -  gI(xj,Pk)uIj'k+1  sgn(nj’fc+1’*)Ax 

3= i 

n 

-  E  Of  («'•'(*„  +M) 

3=1 

n 

~  Y  Dh  Pk’1)  -  g\Xj ,  Pfc))  wj’A'+1]  sgn(uj’fc+1’*)Aa; 


i=i 


<  /’*(0,P^kW1  +  sup  | g^x^P^-g^x^P 


k\ |  ||  7,fc+l|| 

U  ||  BV 


l<j<n 


+  sup  |D^  (s,’‘(^,P*’‘)-9,(iJ,P*))|  (Halloo +  (ll«'*tlMi 

l<ji<n 

By  (2.7),  we  have 

/*(o,p^)Kb’fc+11 

<  |^(0,  Pk '*)  -  g\ 0,  Pk )  |  «o’fc+1  +  1^(4)  -  C\tk 


N 


N 


+cai  Y  \\vJ’k’%  +  max  sup  |/3J’*(^-,  PM)  -  /3J(^-,  Pfc)|  Y  ll«J’fe||i- 
j=i  *<  i<i<™  J=1 

Summing  (2.8)  over  /  =  1,  2, . . . ,  N,  and  using  (2.9)  and  (2.10),  we  obtain 


N 


N 


Z)  IK,fe+M||i  -  E  ||u7’fc,<iii 

i=i  i=i 


At 


N 


(2.9) 


(2.10) 


<  max  sup  |ZV(9«(ij,P*l)-/(xj,.P‘))|  IV  max  J  .V‘+1U  +  £  ||mj  „ 


,hfc+l| 


1<I<N  i <j<n 


l(xj,  Pk'1)  —  g\Xj,  Pk)  max  llz/7,  \\Bv 

J  J  1<I<N 


1<I<N 
,I,k+ 1 1 


1=1 


+N  max  sup  J 

1<I<N  i <j<„  1 

+1V  max  |g7,l(0,  Pfe,?)  —  </(0,  Pk)\  max  ||z/7’fc+1||  yv  max  I CI,l(tk)  —  C1  (tk 

1<I<N  '  '  I  11  "  I^T^AT  I  v  v 


1<I<N 


1<I<N 


N 


N 


+N  maxN  sup  |/3J’*(^-,  Pfc’*)  -  /3J(^-,  Pfc)  |  Y  +  Nuji  J2 


j= i 
v 


7=1 


+  max  sup  |m7’*(xj,  Pk,t)  —  mI(xj,  Pk)\ 
1<I<N  l<j<n 

Noticing  that 


7=1 


|9'-1(i,,Pt'1)-g,(x„Ft)| 

<  l/'Hi,,^)  -9«(i,,Pl)|  +  -/(*,, P*)|, 


12 


we  have  from  (H4)  the  following: 


Similarly,  we  can  show  that 


and 


max  sup  I /3I’t(xj,Pk’'t)  —  f3I(xj,Pk) 

1<I<N  l<j<n 


<  Li 


N 

£ 

7=1 


||n 


I,k,i  I 


i  +  max 

KI<N 


sup  W^Xj^pb)  -  (5\Xj,Pk) 

l<j<n 


max  sup  Pk,i )  —  m 1  ( Xj ,  Pk ) 


1  <I<N  i <j<n 
N 


< 


L2J2\\vI,k,%+  max  sup  |  Pk)  -  m1  Pk)  . 


i=i 


1  <I<N  i <j<n 


Furthermore,  straightforward  computations  yield 


\D~h  [g,-'(xj,Pk--)-g'(x„Pt)]\ 

hit  y^i+H-r^-uP^Ar 


-1  d 


^g^rxj  +  (1  -  r)xj-i,  Pk 


9x\rxj  +  (1  —  r)xj- 1,  P  ,l)dr  —  /  gx{rxj  +  (1  —  r)xj-±,  P  )dr 


< 


g^irxj  +  (1  -  r)xj- 1,  Pk,i )  -  g^irxj  +  (1  -  r)xj- 1,  Pk )  dr 


+  /  +  (1  -  P  )  -  gx(rxj  +  (1  -  r)xj_i,  P  )|  dr. 

J  o 

Hence,  from  (H4)  we  obtain 

max  sup  | D~  [g^xj,  Pk’1)  -  g\xj,  Pfc)]  | 


1<I<N  i <j<n 
N 


< 


Li  H^ll1  +  m^  sup  /  \g^l{xj,  Pk )  -  giixj,  Pk )  |  dr , 


7=1 


1<I<N  i <j<n  J o 


where  ay-  =  ray-  +  (1  —  r)ar,_i.  Set 


iv 


AT 


4  =  L3  N  max  Halloo +  V||u— 111 

-\<T<N  1  1  ^ ' 


,I,k+ 1 1 


7=1 


7=1 


+iVL3  f  max  ||7/’fc+1||w  +  max  ||m/’A:+1||0o1  +£2^ 

\1<I<N  1<T<N  / 


1<7<1V 


|a/’fc+1||1 


7=1 


and 


Pk,i 


N  max  ||  u 

1<I<N 


I,k+ 1 1 


N 

E 

7=1 


U 


I,k+ 1 1 


max  sup  [  | Pk )  -  g*x(xj,  Pk )  |  dr 

1</<7V  i<j<„  Jo 


+1V  max  ||M7’fc+1||_Bvr  max  sup  I gI,l(xj,  Pk)  —  g1  (xj,  Pk 

1</<7V  l<j<n 


+N  max  ||«7’fc+1||oo  max 

1<I<N  1<I<N 


gI,i(o,Pk) 


g\^Pk 


+  N  max 

1<I<N 


c\tk ) 


N 


+NYJ\wLk 

i=i 


i  max 

1<I<N 


sup  \pI’i(xj,Pk)-pI(xj,Pk) 

l<j<n 


N 


+  llw/’fc+1  II1  max  sup  I Pk )  -  m^Xj,  Pk) 


i=i 

Then,  we  have 


1<I<N  i <j<n 


N 

J2  ||u7’fc+M| 

1=1 


N 

i  —  Y2  ||u/’fc,?  | 
i=i 


At 


N 

—  $k  ^  ^  Ill’ll  +  Pk,i- 
1=1 


Since  for  each  k,  pk,i  — >  0  as  i  — >  oo,  the  desired  result  easily  follows  from  this  inequality. 


Theorem  2.4  Suppose  that  Qm  is  a  sequence  of  compact  subsets  of  Q.  Moreover,  assume 
that  for  each  q  e  Q,  there  exists  a  sequence  of  qM  £  Qm  such  that  qM  —>  q  as  M  —>  oo.  Then 
the  functional  SfAx.At  has  a  minimizer  over  Qm-  Furthermore,  if  qlM  denotes  a  minimizer  of 
ST  Ax,.  At,  over  Qm  and  Axi,  At*  — >  0,  then  any  subsequence  of  qlM  has  a  further  subsequence 
which  converges  to  a  minimizer  of  J . 

Proof.  The  proof  of  this  theorem  is  a  direct  application  of  the  abstract  theory  in  [18],  based 
on  the  convergence  of  — ►  J{q)-  ■ 


3  Numerical  Results 

In  this  section,  we  present  ample  numerical  simulations  and  statistical  results.  In  all  of  the 
simulations  below  we  assume  L  —  1,  T  =  1,  and  CT(t)  =  0  for  /  =  1,  2, . . . ,  N. 

In  subsections  3.1  and  3.2,  we  assume  N  =  1  and  that  all  the  parameters  are  known 
except  for  f3.  To  estimate  f3  we  use  data  which  are  generated  computationally  as  follows: 
Let 


u°(x)  =  3exp(— 2(x  —  0.5)2),  g(x,  P )  =  5(1  —  x)  exp(— 3P), 

m{x ,  P)  =  exp(4(a:  —  0.4)2)  exp(0.2P),  fi(x,  P)  =  6x(l  —  x)  exp(— 3P), 
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and  we  solve  (2.1)  and  (2.3)  for  UAXiAt(x,t)-  We  set  the  data  Zk  =  (l+efc)  J  UAx,At(xAk)dx, 
where  ek  is  a  random  sample  from  a  normal  random  number  generator  with  mean  zero  and 
standard  deviation  a  =  0.02. 


3.1  1  —  D  linear  estimation  problem  for  finite  dimensional  param¬ 

eter  space  when  N  =  1 

In  our  first  example  we  assume  that  (3  is  of  a  separable  form  given  by  f3(x,  P)  =  b(x)  exp(— 3 P), 
where  b(x)  =  fj,x(l  —  xv )  with  /i  and  v  two  unknown  constants  to  be  identified.  Hence,  the 
solution  to  our  least-squares  problems  involves  identifying  the  two  constants  /i  and  v  from  a 
compact  subset  of  so  as  to  minimize  the  least-squares  cost  functional 


M)  = 


m  ,  r\ 

E  >°s  (/ 

4.1  VJ0 


UAx,At(x ,  tk)  (-l)dx  +  11-  log (Zk  +  1) 


In  order  to  test  the  performance  of  the  parameter-estimation  technique  when  no  infinite 
dimensional  effects  are  present,  in  Figure  1  we  choose  Ax  =  At  =  0.005  for  both  generating 
the  data  and  the  numerical  solution  (2.3)  in  the  least-squares  problem.  This  avoids  the 
infinite-dimensional  effect  of  the  partial  differential  equation  given  in  (1.1).  In  fact,  if  the 
noise  is  removed  from  the  data,  and  the  parameters  /i  and  v  are  known,  then  numerically 
solving  our  model  produces  the  exact  data. 

In  Figure  2  we  use  Ax  =  At  =  0.005  to  generate  the  data  while  we  use  Ax  =  At  =  0.01 
for  the  numerical  solution  (2.3)  in  the  least-squares  problem.  Thus,  in  this  case  the  data  are 
not  exactly  attained  by  our  model  even  if  the  noise  is  removed  (an  error  is  present  due  to  the 
finite-dimensional  approximation  of  our  infinite-dimensional  model).  The  results  of  Figure 
2  are  obtained  by  using  the  same  values  for  the  rest  of  the  parameters  as  those  of  Figure  1. 


A  similar  format  for  presenting  the  results  of  1000  inverse  problem  calculations  was  used 
in  Figure  1  and  2.  The  left  part  of  each  of  the  figures  represents  the  S  (for  our  case  S  =  1000) 
numerical  results  for  the  estimated  parameter  bs  (x)  ( s  =  1,  2, . . . ,  S)  versus  the  exact  b(x), 
where  these  1000  distinct  numerical  results  graphed  were  obtained  by  solving  1000  inverse 
problems,  each  of  which  corresponds  to  a  given  noise  sample  {e*,}.  The  right  part  represents 
the  figure  of  the  corresponding  95%  confidence  interval  (dashed  line)  versus  the  exact  b(x) 
(solid  line),  where  the  95%  confidence  interval  is  obtained  by  choosing  the  band  between 


the  upper  2.5%  and  lower  2.5%  of  these  1000  numerical  results.  Table  1  provides  statistical 

1S 

results  for  the  corresponding  graphs,  where  AB(x)  =  —  Ub‘  (x)  —  b(x))  denotes  the  average 

S= 1 

AB(x) 

bias  for  all  approximations  at  x,  RAB(x )  =  100 — - denotes  the  relative  average  bias  for 

b(x) 
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all  approximations  at  x  and 


SE(x) 


J2(bs(x)  -  b(x)  -  AB(x))2 

S= 1 


1 

S-l 


s 


X>*w 


denotes  the  sampling  standard  error  at  the  point  x.  Note  that  this  is  simply  the  usual 
asymptotic  formula  for  the  pointwise  standard  error  (e.g.,  see  p.  28,  37  of  [21]  and  p.  308  of 
[45]). 

Although  the  estimates  in  both  figures  are  good,  the  results  in  Figures  1-2  and  Table 
1  suggest  that  infinite- dimensional  effects  can  lead  to  a  slightly  under  biased  estimator. 
We  suspect  that  this  bias  depends  on  the  choice  of  the  numerical  scheme  used  for  solving 
the  infinite-dimensional  partial  differential  equation  model.  Here  we  are  using  an  upwind 
scheme  for  approximating  the  model  and  a  right-hand  sum  for  approximating  all  the  integrals 
involved.  This  biased  estimator  may  be  improved  if,  for  example,  a  centered  finite  difference 
approximation  is  used  together  with  a  trapezoidal  rule  for  integration. 


Figure  1:  Ax  =  At  =  0.005  to  generate  the  data  and  solve  the  least-squares.  For  the  left 
part  of  the  figure,  each  of  the  grey  lines  (....)  denotes  a  distinct  result  for  a  given  sample 

{**}■ 


The  above  statistical  results  (essentially  on  how  measurement  error  affects  estimates)  are 
based  on  a  large  number  of  numerical  simulations  (somewhat  in  the  spirit  of  Bayesian  based 
MCMC  calculations  used  to  estimate  means  and  variances  in  a  probability  distribution  from 
“experimental”  data).  Any  estimate  of  model  parameters  from  data  can  also  be  accompanied 
by  an  estimate  of  uncertainty  using  standard  regression  formulations  from  statistics  [21]. 
Thus,  in  the  remaining  part  of  this  subsection,  we  present  a  statistical  based  method  to 
actually  compute  the  variance  in  the  estimated  model  parameters  q  =  (/i,  v). 
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Figure  2:  Ax  =  At  =  0.005  to  generate  the  data  and  Ax  =  At  =  0.01  to  solve  the  least- 
squares.  For  the  left  part  of  the  figure,  each  of  the  grey  lines  (....)  denotes  a  distinct  result 
for  a  given  sample  {e/c}. 


X 

AB(x) 

RAB(x) 

SE(x) 

X 

AB(x) 

RAB(x) 

SE(x) 

0.1 

-0.0037 

-0.6870 

0.0749 

0.1 

-0.0390 

-7.2314 

0.0747 

0.2 

-0.0092 

-0.9580 

0.0993 

0.2 

-0.0651 

-6.7812 

0.1053 

0.3 

-0.0107 

-0.8463 

0.0975 

0.3 

-0.0768 

-6.0949 

0.1130 

0.4 

-0.0079 

-0.5497 

0.0860 

0.4 

-0.0763 

-5.2995 

0.1124 

0.5 

-0.0021 

-0.1427 

0.0798 

0.5 

-0.0666 

-4.4422 

0.1138 

0.6 

0.0049 

0.3378 

0.0852 

0.6 

-0.0511 

-3.5460 

0.1188 

0.7 

0.0110 

0.8707 

0.0926 

0.7 

-0.0331 

-2.6236 

0.1202 

0.8 

0.0138 

1.4425 

0.0882 

0.8 

-0.0162 

-1.6830 

0.1075 

0.9 

0.0110 

2.0444 

0.0605 

0.9 

-0.0039 

-0.7294 

0.0706 

Table  1:  Left  and  right  tables  are  statistical  results  for  Figure  1  and  Figure  2,  respectively. 
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To  perform  this  analysis,  we  need  to  compute  the  sensitivity  matrix 


X(q) 


PM(U;g) 

Pu(ti,q) 

1  +  P(U;  q) 
Pfifaq) 

1  +  P(U;  q) 
Pvfaq) 

1  +  P(t2;  q) 

1  +  P{t2 ;  q) 

Pfj.it mi  q) 

Pyitmi  q) 

1  +  P(Urh  q) 

l  +  Pitm'i  q) 

(3.1) 


Note  that  we  cannot  compute  P(t;q),  P^(t]q)  and  Pv{t\  q)  directly  from  our  model.  There¬ 
fore,  we  use  the  difference  scheme  (2.1)  to  obtain  the  following  approximation  of  P(t ;  q): 


PM 


UAx,At(x,ti(l)dx. 


Then  we  use  a  forward  difference  approximation  for  the  derivative  P/i(f;  q)  and  Pu(t;  q)  given 
by 

P^(t;  /x,  v)  =  —  /x  +  A/x,  v)  -  P{t ;  /x,  z/)) 

and 

Pu{t]  q)  =  ^  (p(U  /U  i'  +  Az/)  -  P(f;  /X,  I/))  . 

Substituting  P{ti,q),  P^q)  and  Pv{ti,q)  for  P(tpg),  P^t^q)  and  Pu{U;q)  in  (3.1),  re¬ 
spectively,  we  obtain  the  following  approximation  of  X{q): 


Pn(ti;q) 

Pu(ti',q) 

i  +  P(U;  9) 
Ppfaq) 

l  +  P(U;  9) 
Pvfaq) 

xiq)  = 

1  +  P(t2;  9) 

1  +  P(t2;  9) 

P flit  mi  q) 

1  +  Pitm'i  q) 

1  +  Pitm'i  q) 

Under  standard  assumptions  of  classical  nonlinear  regression  theory,  we  know  that  if 
Q  ~  A/"(0,a2),  where  e,-  is  the  difference  between  observation  and  model  at  time  tt,  then 
the  least-squares  estimate  q*  is  expected  to  be  asymptotically  normally  distributed.  In 
particular,  for  large  samples,  we  may  assume 


9 *  ~  Af[g0,  a2{XTiq0)X{q0)}-1],  (3.2) 

where  q0  is  the  true  vector  of  parameters  and  a2{A"T(g0)AT(g0)}-1  is  the  true  covariance 
matrix  (see  [21],  Chapter  2). 
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Since  go  and  a2  are  not  available,  we  follow  a  standard  statistical  practice  [5] :  substitute 
the  computed  estimate  g*  for  g0  and  approximate  a2  by 

(log  (-P(0;9*)  + 1)  -  log (z,  + i>)  (3.3) 

3= 1 

in  (3.2)  to  obtain  the  standard  deviation  for  our  estimates.  In  particular,  if 

v'  =  =  [  l"11  l)2  ] , 

r  21  V22 

then  we  take  a/^ii  and  -\/V^2  to  be  the  standard  deviation  for  parameters  /i  and  z/,  respec¬ 
tively.  The  following  two  tables  are  the  standard  deviation  of  /i  and  v  for  the  results  of  the 
first  eight  numerical  simulations  of  Figure  1  and  Figure  2,  respectively. 


h 

1.1613 

1.0494 

1.0451 

1.1109 

1.0864 

1.4684 

1.1605 

1.0512 

V 

1.2124 

0.3073 

0.2999 

0.2741 

0.2701 

1.5555 

0.2482 

0.2390 

Table  2:  Standard  deviation  for  the  results  of  the  first  8  numerical  simulations  of  Figure  1. 


l< 

1.7066 

1.5636 

1.6192 

1.7974 

1.6389 

2.8009 

1.8619 

1.3893 

V 

0.7716 

0.3238 

0.4838 

0.1812 

0.3426 

2.8685 

0.3828 

0.4136 

Table  3:  Standard  deviation  for  the  results  of  the  first  8  numerical  simulations  of  Figure  2. 

Table  4  provides  the  average  standard  deviation  of  (i  and  v  for  the  results  of  all  the  1000 
numerical  simulations  of  Figure  1  and  Figure  2,  respectively.  We  note  that  in  most  practical 
situations  using  experimental  data,  one  does  not  expect  to  have  1000  experiments  performed. 
But  the  above  procedures  will  produce  estimates  of  variances  even  in  the  case  when  one  has 
only  one  data  set! 


Figure  1 

Figure  2 

h 

1.1921 

1.9197 

V 

0.4566 

0.8572 

Table  4:  Average  of  standard  deviation  for  all  the  results  of  the  numerical  simulations  of 
Figures  1-2. 

3.2  1  —  D  linear  estimation  problem  for  infinite  dimensional  pa¬ 

rameter  space  when  N  =  1 

In  this  example,  we  assume  that  (3  is  of  a  separable  form  given  by  /3(x,  P)  =  b(x)  exp(— 3 P), 
where  b(x)  is  an  unknown  parameter  that  we  want  to  identify. 
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Let 


v  =  {fe  C[0, 1]  :  I  f(x)  -  f(y) I  <  K\x  -  y\,f(0)  =  /( 1)  =  0}. 

Choose  the  parameter  space  Q  =  T>.  Clearly,  by  Arzela-Ascoli  Theorem  [33]  Q  is  compact 
in  C[0, 1].  We  approximate  the  infinite  dimensional  parameter  space  as  follows:  For  M  a 
positive  integer  and  b  £  Q,  we  set 

M— 1 

(IMb)(x)  =  b 

i= 1 


where  (f)lM(x,0, 1)  are  the  linear  spline  functions  on  a  uniform  mesh  of  the  interval  [0,1]. 
These  are  defined  by 


<Pm{x]  0,1)  =  < 


.  X 

1  -  1  +  7, 

h 

(i  —  1  )h  <  x  <  ih, 

.  •  x 

ih  <  x  <  (i  +  1  )h, 

i  =  1,  2, . . . ,  M  —  1, 

0, 

\x  —  ih  >  h, 

where  h  =  — .  It  can  be  readily  argued  that  lim  Xa/6  =  6  in  C[0, 1],  uniformly  in  6  [44], 
M  M— >oo 

Hence,  if  bM  G  Qm  =  %m(Q)  is  given  by 


M— 1 

Mar)  =  Am0m(x;M), 

i=l 


then  the  solution  of  our  finite  dimensional  identification  problem  involves  identifying  the 
M  —  1  coefficients  &'om  a  compact  subset  of  M^_1  so  as  to  minimize  the  least- 

squares  cost  functional  (2.4). 

In  order  to  indirectly  implement  the  compactness  constraints  of  Q ,  we  use  a  regularized 
least  squares  cost  functional  of  the  form 


J\x.aM)  =  V]  los  (  / 

k= 1  A 


UAx,At(x,  4;  q)dx  +  11-  log (Zk  +  1) 


2  f1 

d 

-rbM[x) 

Jo 

dx  v  ; 

dx, 


where  a  >  0  is  the  regularization  parameter. 

The  left  part  of  each  of  the  following  figures  again  represents  the  S  (=1000)  numerical 
results  of  the  estimated  parameter  versus  the  exact  parameter  b(x).  The  right  part  repre¬ 
sents  the  figure  of  the  corresponding  95%  confidence  interval  (dashed  line)  versus  the  exact 
b(x)  (solid  line).  The  tables  provide  statistical  results  for  the  corresponding  graphs. 
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Effect  of  infinite- dimensional  model  on  parameter  estimate. 


In  Figure  3  we  use  Aa;  =  0.005  and  At  =  0.005  to  generate  the  data  and  the  numerical 
solution  (2.3)  for  the  least-squares  problem.  This  removes  the  infinite-dimensional  effect  of 
the  partial  differential  equation  given  by  (1.1).  However,  in  Figure  4  we  use  Ax  =  At  =  0.005 
to  generate  the  data  and  Ax  =  At  =  0.01  to  compute  (2.3).  Thus,  in  this  case  the  data 
are  not  exactly  attained  by  our  model  even  if  the  noise  is  removed.  We  observe  that  while 
the  estimates  in  both  figures  are  good,  the  results  in  Figures  3-4  and  Table  5  suggest  that 
infinite-dimensional  effects  can  lead  to  a  slightly  under  biased  estimator. 


Figure  3:  M  =  10,  a  =  3e  —  5.  Each  of  the  grey  lines  (....)  of  the  left  part  of  the  figure 
denotes  a  distinct  result  for  a  given  sample  {e*;,}. 


X 

AB(x) 

RAB(x) 

SE(x) 

X 

AB{x) 

RAB(x) 

SE(x) 

0.1 

-0.0778 

-14.4108 

0.0723 

0.1 

-0.1236 

-22.8940 

0.0667 

0.2 

-0.0816 

-8.5015 

0.1070 

0.2 

-0.1571 

-16.3628 

0.1040 

0.3 

-0.0400 

-3.1727 

0.1012 

0.3 

-0.1284 

-10.1885 

0.1141 

0.4 

0.0110 

0.7636 

0.0834 

0.4 

-0.0785 

-5.4485 

0.1130 

0.5 

0.0386 

2.5745 

0.0818 

0.5 

-0.0440 

-2.9329 

0.1110 

0.6 

0.0283 

1.9621 

0.0868 

0.6 

-0.0446 

-3.0966 

0.1049 

0.7 

-0.0124 

-0.9880 

0.0779 

0.7 

-0.0754 

-5.9875 

0.0885 

0.8 

-0.0559 

-5.8206 

0.0556 

0.8 

-0.1059 

-11.0334 

0.0624 

0.9 

-0.0623 

-11.5426 

0.0280 

0.9 

-0.0939 

-17.3949 

0.0323 

Table  5:  Left  and  right  tables  are  statistical  results  for  Figure  3  and  Figure  4,  respectively. 


Effect  of  regularization  parameter  a  on  parameter  estimate. 
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Figure  4:  M  =  10,  a  =  3e  —  5.  Each  of  the  grey  lines  (....)  of  the  left  part  of  the  figure 
denotes  a  distinct  result  for  a  given  sample  {e^}. 


In  Figures  5  and  6  we  change  the  parameter  a  while  keeping  the  rest  fixed.  Clearly,  low 
regularization  parameter  leads  to  relatively  bad  estimates  although  the  estimator  in  this  case 
seems  to  be  the  least  biased  (see  Figure  5  and  left  part  of  Table  6).  Increasing  the  value 
of  a  leads  to  better  parameter  estimates,  but  the  estimator  becomes  more  under  biased 
(see  Figure  6  and  right  part  of  Table  6).  If  this  value  is  increased  more,  the  estimator  is 
more  biased.  Also  the  parameter  estimate  becomes  worse  than  before.  This  suggests,  not 
surprisingly,  that  there  is  an  optimal  choice  for  the  parameter  a  which  produces  the  best 
results  for  the  parameter  estimates. 


Figure  5:  M  =  10,  a  =  le  —  5.  Each  of  the  grey  lines  (....)  of  the  left  part  of  the  figure 
denotes  a  distinct  result  for  a  given  sample  {e/c}. 
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X 

AB(x) 

RAB(x) 

SE(x) 

X 

AB(x) 

RAB(x) 

SE(x) 

0.1 

-0.1277 

-23.6389 

0.1206 

0.1 

-0.1241 

-22.9816 

0.0506 

0.2 

-0.1648 

-17.1644 

0.1791 

0.2 

-0.1621 

-16.8881 

0.0842 

0.3 

-0.1284 

-10.1938 

0.1618 

0.3 

-0.1432 

-11.3627 

0.1011 

0.4 

-0.0599 

-4.1591 

0.1221 

0.4 

-0.1050 

-7.2906 

0.1078 

0.5 

-0.0072 

-0.4806 

0.1169 

0.5 

-0.0791 

-5.2736 

0.1087 

0.6 

0.0026 

0.1788 

0.1274 

0.6 

-0.0837 

-5.8139 

0.1009 

0.7 

-0.0253 

-2.0101 

0.1126 

0.7 

-0.1077 

-8.5443 

0.0847 

0.8 

-0.0631 

-6.5678 

0.0780 

0.8 

-0.1288 

-13.4165 

0.0602 

0.9 

-0.0642 

-11.8944 

0.0427 

0.9 

-0.1042 

-19.3027 

0.0313 

Table  6:  Left  and  right  tables  are  statistical  results  for  Figure  5  and  Figure  6,  respectively. 


Figure  6:  M  =  10,  a  =  5e  —  5.  Each  of  the  grey  lines  (....)  of  the  left  part  of  the  figure 
denotes  a  distinct  result  for  a  given  sample  {efc}. 
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3.3  1  —  D  linear  estimation  problem  for  infinite  dimensional  pa¬ 

rameter  space  when  N  =  2 


In  this  section,  we  assume  N  —  2  and  that  all  the  parameters  are  known  except  for  (3 1 
and  (32 .  To  estimate  (3 1  and  /32,  we  assume  that  they  are  of  a  separable  form  given  by 
(31{x1P)  =  61(x)  exp(— P)  and  /32(x,P )  =  b2(x)  exp(— P),  respectively,  where  61(:r)  and 
b2(x )  are  unknown  parameters  to  be  identified.  To  estimate  bl(x)  and  b2(x),  we  use  data 


which  are  generated  computationally  as  follows:  Let  r)I'J  = 


1, 

0, 


I  =  J 


for  Figure  7  and 


ryIJ  =  0.5,  I,J  =  1,2  for  Figure  8,  uI,0(x )  =  3exp(— 2(x  —  0.1)2),  and  for  the  parameters 
g1  ,  'm1  and  (3 1  we  use  the  following  choice  of  functions: 


g 1  =  2(1  —  x)  exp(— 0.8P),  g2  =  (1  —  rr)(l  +  2 P)  exp (— P), 

m1  =  exp(2(a:  —  0.4)2)  exp(0.2P),  m2  =  exp(2(a:  —  0.4)2)  exp(0.2P), 

/ 3 1  =  6(1  —  x)a;exp(— P),  (3 2  =  6(1  —  rr)a;exp(— 5(x  —  0.5)2)  exp(— P), 


and  solve  (1.1)  for  U^x  At(x,t),  I  —  1,2.  We  set  the  data  =  (1  +£itk  )  uL,a  t(x,tk)dx, 


2  /»1 

7  =  1,2  for  Figure  7  and  Z*,  =  (1  +  £k)  /  UAx  At{x,  tk)dx  for  Figure  8,  where  and 

/=i  do 

both  are  the  random  sample  from  a  normal  random  number  generator  with  mean  zero  and 
standard  deviation  a  =  0.02. 

We  choose  the  parameter  space  Q  =  D  x  V.  Clearly,  Q  is  compact  in  C[ 0, 1]  x  C[0,1], 
We  approximate  the  infinite  dimensional  parameter  space  as  follows:  For  M1;  M2  positive 
integers  and  any  (61,62)  £  Q,  we  set 


(lMjbJ)(x) 


Mj—1  /  .  \ 

E ',j(P)«jj(:c;0'1)'  J  =  1-z 

i=  1  ' 


Clearly,  lim  lMj b'J  =  bJ  in  C [0, 1] ,  uniformly  in  6J,  J  =  1,2.  Hence,  if  bJM  e  Qmj  = 

ATj— >00  J 

JMj(<5)  is  given  by 

Mj-l 

bJMj(x)  —  ^Mj4>lMj{x't  0, 1),  J  =  l,2, 

then  the  solution  of  our  finite  dimensional  identification  problem  involves  identifying  the 
Mi  +  M2  —  2  coefficients  {X Mj}i=i~j=i  from  a  compact  subset  of  ]g^i+M2-2  so  as  to  minimize 
the  least-squares  cost  functional  (2.4)  or  (2.5). 

In  order  to  indirectly  implement  the  compactness  constraints  of  Q,  we  still  use  the  regu- 
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larized  least-squares  cost  functional.  For  Figure  7  we  use  the  form 


2  in  /  i 

jaxMq)  =  5Z5Z  log( 

1=1  k= i  'v"'° 


uL.At{xi  4;  q)dx  +  11-  log(z/ifc  + 1) 


X> 


^bMj(X) 


1=1 


1 0 


dx 


dx, 


and  for  Figure  8  we  use  the  form 

m 

Ax,a t(o)  ^  ^ 

=1 

E 


k= 1 


iog(t/4. 

\i= i 


4;  +  1  -  log(Zfc  +  1) 


2  r 

OLI 


1=1 


'0 


dx, 


where  a/  >  0,  7  =  1,2  are  the  regularization  parameters  and  m  =  100  for  Figures  7  and  8. 

In  the  rest  of  our  simulations  we  use  Ax  =  A t  —  0.005  to  generate  the  data  and  Ax  = 
At  =  0.01  to  solve  the  least-squares.  Thus,  in  these  cases  the  data  are  not  exactly  attained 
by  our  model  even  if  the  noise  is  removed. 

The  upper-left  part  and  the  lower-left  part  of  the  following  two  figures  represent  the  S 
(=1000)  numerical  results  of  the  estimated  parameters  blMl(x)  and  b2M2(x)  versus  the  exact 
parameters  61(a:)  and  b2{x),  respectively.  The  upper- right  part  and  the  lower  right  part 
represent  the  figures  of  the  corresponding  95%  confidence  interval  (dashed  line)  versus  the 
exact  61(x)  and  b2{x)  (solid  line),  respectively.  The  tables  provide  statistical  results  for  the 
corresponding  graphs. 

Note  that  the  results  in  Figure  7  and  Table  7  are  slightly  better  than  those  in  Figure  8 
and  Table  8.  This  is  expected  since  in  Figure  7  we  are  sampling  data  for  each  of  the  two 
populations,  which  provides  more  information  than  sampling  the  sum  of  the  two  populations 
only,  as  is  the  case  in  Figure  8.  Also  note  that  in  both  of  these  figures  we  let  M  =  M\  = 
M2  =  10. 
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Figure  7:  M  =  10,  a.\  =  5e  —  5,  a-2  =  5e  —  5.  Each  of  the  grey  lines  (....)  of  the  left  part  of 
the  figure  denotes  a  distinct  result  for  a  given  sample  {£/,&}• 


X 

AB(x) 

RAB(x) 

SE{x) 

X 

AB{x) 

RAB(x) 

SE(x) 

0.1 

-0.0187 

-3.4717 

0.0880 

0.1 

0.1684 

69.4034 

0.0959 

0.2 

-0.0004 

-0.0447 

0.1276 

0.2 

0.1628 

26.5887 

0.1528 

0.3 

0.0334 

2.6514 

0.1053 

0.3 

0.0487 

4.7244 

0.1483 

0.4 

0.0562 

3.9007 

0.0493 

0.4 

-0.0728 

-5.3114 

0.0946 

0.5 

0.0449 

2.9941 

0.0548 

0.5 

-0.1134 

-7.5604 

0.0464 

0.6 

-0.0040 

-0.2805 

0.0860 

0.6 

-0.0437 

-3.1871 

0.0860 

0.7 

-0.0683 

-5.4239 

0.0836 

0.7 

0.0931 

9.0282 

0.1053 

0.8 

-0.1101 

-11.4644 

0.0576 

0.8 

0.2039 

33.3052 

0.0819 

0.9 

-0.0929 

-17.2091 

0.0272 

0.9 

0.1954 

80.5164 

0.0402 

Table  7:  Left  and  right  tables  are  statistical  results  of  bl(x )  and  b2(x)  for  Figure  7,  respec¬ 
tively. 
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approximate 

exact 


0  0.1  0.2  0.3  0.4  0.5 


0.7  0.8  0.9  1 


Figure  8:  M  =  10,  ct\  =  5e  —  5,  «2  =  5e  —  5.  Each  of  the  grey  lines  (....)  of  the  left  part  of 
the  figure  denotes  a  distinct  result  for  a  given  sample  {e*,}. 


X 

AB{x) 

RAB(x) 

SE{x) 

X 

AB{x) 

RAB(x ) 

SE(x) 

0.1 

-0.0687 

-12.7279 

0.0765 

0.1 

0.1926 

79.3867 

0.1066 

0.2 

-0.0790 

-8.2334 

0.1096 

0.2 

0.2106 

34.4018 

0.1757 

0.3 

-0.0572 

-4.5419 

0.0891 

0.3 

0.1187 

11.5041 

0.1784 

0.4 

-0.0402 

-2.7920 

0.0435 

0.4 

0.0178 

1.2960 

0.1208 

0.5 

-0.0537 

-3.5800 

0.0588 

0.5 

-0.0069 

-0.4598 

0.0549 

0.6 

-0.0980 

-6.8075 

0.0871 

0.6 

0.0665 

4.8565 

0.0915 

0.7 

-0.1490 

-11.8273 

0.0846 

0.7 

0.1889 

18.3112 

0.1157 

0.8 

-0.1694 

-17.6443 

0.0596 

0.8 

0.2704 

44.1765 

0.0915 

0.9 

-0.1255 

-23.2483 

0.0296 

0.9 

0.2239 

92.2680 

0.0459 

Table  8:  Left  and  right  tables  are  statistical  results  of  bl(x )  and  b2(x)  for  Figure  8,  respec¬ 
tively. 
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4  Concluding  Remarks 


In  this  paper  we  have  developed  a  numerical  technique  for  identifying  unknown  parameters 
in  a  general  size-structured  population  model.  A  main  focus  of  the  paper  is  on  a  statistical 
study  of  the  parameter  estimation  technique.  This  was  carried  out  by  calculating  pointwise 
standard  errors  on  the  estimated  parameters  (functions)  via  use  of  thousands  of  numerical 
experiments. 

Several  conclusions  can  be  drawn  from  our  studies.  1)  The  method  discussed  above  seems 
to  perform  well  and  produce  good  confidence  intervals  for  the  parameters.  2)  When  the 
infinite  dimensional  effects  of  the  model  and  the  parameter  space  are  removed,  the  resulting 
numerical  and  statistical  values  suggest  that  the  least-squares  technique  produces  very  good 
unbiased  parameter  estimates.  3)  The  type  of  numerical  scheme  used  for  approximating 
the  infinite-dimensional  model  as  well  as  the  parameter  space  may  influence  the  bias  in  the 
parameter  estimation  technique.  4)  The  commonly  used  regularization  term  is  crucial  for 
enforcing  compactness  and  obtaining  better  estimates.  However,  it  may  also  introduce  more 
bias  in  the  estimator. 

We  note  in  closing  that  the  system  (1.1)  investigated  in  this  paper  is  a  special  case 
of  the  measure  dependent  aggregate  dynamics  problems  formulated  in  [6]  wherein  individ¬ 
ual  (uncoupled)  dynamics  are  not  available.  Inverse  problems  for  such  systems  have  been 
investigated  in  a  number  of  applications  including  cellular  level  HIV  modelling  [7],  hys¬ 
teresis  in  viscoelastic  materials  [8,  9],  shear  waves  in  biotissue  [10],  and  electromagnetic 
interrogation  in  complex  materials  [11].  In  a  more  general  formulation  (currently  under 
investigation  by  the  authors),  one  has  a  probability  distribution  F  of  individual  parame¬ 
ters  q(x,P )  =  q  =  (g,m,  /3,C)  on  an  admissible  set  Q.  The  system  (1.1)  is  replaced  by  a 
continuum  of  systems  for  u(x,t;q(x,  P))  with  the  total  population  P(t;F )  given  by 


P(f;F) 


/ 

-  r-L 

/  u(x,t',q)dx 

dF(q)  =  I 

-  r-L 

/  u{x,t-,q)dx 

>Q 

Jo 

JQ 

Jo 

the  latter  equality  holding  if  F  has  a  density  /.  The  aggregate  dynamics  for  u  depend 
explicitly  on  F  through  the  dependence  of  the  individual  rate  parameters  (g,  m,  (3,  C )  on  the 
total  population  P. 

If  F  is  a  discrete  measure  with  N  atoms  at  qJ  of  mass  fj ,  then  we  have 


N  fL 


Pit i  F)  =  ±h[ 

J=1 


u(x,  t;  qJ)dx. 


Moreover,  if  F  is  uniformly  and  discretely  distributed  (fj 


this  becomes 


Pit-  F ) 


1 

N 


N  i>L 

/  u(x,t-,qJ)dx, 
j=i  Jo 


which  is  simply  a  scaled  (by  — )  version  of  (1.2).  Of  course,  even  in  this  simple  case,  the 
system  does  not  decouple,  (i.e.,  individual  dynamics  are  not  available).  This  will  be  the  case 
anytime  the  individual  parameters  for  subpopulations  depend  on  the  total  population.  It  is 
also  clear  that  inverse  problems  with  such  measure  dependent  dynamics  are  a  generalized 
version  of  the  estimation  problems  discussed  in  the  statistical  literature  in  the  context  of 
hierarchial  or  mixed  effects  modelling  [20,  21,  22], 
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